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Abstract. Dynamics of the repulsive Bose-Einstein condensate (BEC) in a double- 
well trap is explored within the 3D time-dependent Gross-Pitaevskii equation. The 
model avoids numerous common approximations (two-mode treatment, time-space 
factorization, fixed values of the chemical potential and barrier penetrability, etc) and 
thus provides a realistic description of BEC dynamics, including both weak-coupling 
(sub-barrier) and strong-coupling (above-barrier) regimes and their crossover. The 
strong coupling regime is achieved by increasing the number N of BEC atoms and 
thus the chemical potential. The evolution with N of Josephson oscillations (JO) 
and Macroscopic Quantum Self- Trapping (MQST) is examined and the crucial impact 
of the BEC interaction is demonstrated. At weak coupling, the calculations well 
reproduce the JO/MQST experimental data. At strong coupling, with a significant 
overlap of the left and right BECs, we observe a remarkable persistence of the 
Josephson-like dynamics: the JO and MQST converge to a high-frequency JO-like 
mode where both population imbalance and phase difference oscillate around the zero 
averages. The results open new avenues for BEC interferometry. 
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1. Introduction 

The trapped Bose-Einstein condensate (BEC) is nowadays widely recognized as a source 
of new fascinating physics, see monographs [TJ |2] as well as early [3J, HI E] and recent 
[6j [7] reviews. Among diverse aspects of this field, a large attention is paid to dynamics 
of bound condensates and relevant nonlinear effects caused by the interaction between 
BEC atoms [8]. A boson Josephson junction in a double- well trap represents a 
typical example of relevant system. Its basic dynamical regimes, Josephson oscillations 
(JO) and Macroscopic Quantum Self- Trapping (MQST), were widely investigated in the 
weak coupling limit, both in theory [El EJ El El] and experiment [HJ [12], [13] . Last years, 
new processes in bound condensates, like controlled BEC transport [HI [To] , matter 
wave interferometry |16j . squeezing and entanglement [17], were in the focus. Being 
diverse and sophisticated, these processes concern, nevertheless, the basic JO and MQST 
dynamics and, what is important, often occur beyond the weak coupling limit. Thus, the 
study of JO/MQST at strong coupling (SC), e.g. in the crossover between sub-barrier 
and above-barrier transfer, becomes essential. 

The SC regime routinely arises at a strong interaction or/and large number of BEC 
atoms. Hence this regime is quite common. At the same time, its features have not been 
yet properly investigated. Note that at SC, the left and right BECs in a double- well trap 
are not well separated, so that the use of the ordinary variables of Josephson dynamics, 
population imbalance z and phase difference 6, is not well justified. There thus arises a 
general question whether the JO and MQST survive at SC and, if yes, to what extent 
they are modified? 

Numerous previous studies of Josephson dynamics in a double-well trap were 
performed within the two-mode approximation (TMA) [18], where only two lowest 
energy levels of the system contribute to the dynamics and a weak coupling of the 
condensates is assumed. Obviously, SC dynamics can involve many energy levels and the 
TMA is then not correct. Even the TMA modifications, like a variable tunneling model 
[19] and a direct implementation of Wannier states |20j, do not amend the principle 
TMA limitation to deal only with two lowest BEC states and thus leave TMA ineligible 
to the SC case. 

The SC dynamics requires a more involved theory embracing impact of all excited 
states and avoiding, as much as possible, other standard approximations, e.g. the space- 
time factorization and partition (for every well) of the total order parameter. There are 
already some studies of this kind [HI [H [21] [221 123] . Most of them consider transformation 
of the trap from a double to a single well shape [21] or back [221 [23] and thus partly 
concern the SC case. The treatment varies from using the Gross-Pitaevskii equation 
[21] in [21] to many-body quantum dynamics in [221 123] . These studies confirm the 
important role played by the excited states, beyond the TMA. However, these studies 
do not especially address the evolution of the Josephson dynamics when approaching 
the SC and do not provide a relevant analysis. Moreover, these studies are usually 
performed in the one-dimensional (ID) limit [25], the relevance of which is not well 



Strong-coupling dynamics of Bose-Einstein condensate in a double-well trap 



3 



demonstrated for the SC. In addition, the ID limit is obviously questionable in a true 
three-dimensional (3D) case, which is often met in practice, see e.g. the JO/MQST 
Heidelberg experiment [TT] [12] . 

In this paper, we analyze the evolution of the JO/MQST dynamics in a double- well 
trap while transforming the system from the weak coupling to the SC. The transfer is 
achieved by increasing the number of BEC atoms from iV=1000 to 10000. We thus 
obtain the rise of the cumulative effect of the repulsive interaction and the subsequent 
upshift of the chemical potential. The analysis is based on an explicit solution of the 
3D time-dependent Gross-Pitaevskii equation for the total order parameter. None of 
the questionable approximations mentioned above is thus used. Both the population 
imbalance z and phase difference 9 are examined. 

The calculations are performed for a double- well configuration and the JO/MQST 
initial conditions of the Heidelberg experiment [HJ IT2] . For iV=1000, we reproduce 
the JO/MQST experimental data [TTJ [12] for the repulsive condensate of 87 Rb atoms. 
This justifies the relevance of our model in the weak coupling limit. By increasing the 
number of atoms up to iV=10000, we approach the SC regime and show that, despite an 
essential intersection of the left and right BECs in the double-well trap, these BECs still 
keep their individuality. Hence, using the relative Josephson conjugate variables z and 
9 remains reasonable. We show that JO survive at SC, though with a higher frequency. 
The MQST is transformed to a similar JO. Actually, JO and MQST merge to the same 
mode. Perhaps, a similar MQST — > JO transfer was earlier predicted for ID BEC as a 
reappearance of tunneling in the strong interaction limit [26J. 

In most of the previous calculations (see, e.g. [5]), the initial conditions are obtained 
by constraining the system into a non-stationary state of the symmetric trap. This may 
be questionable at SC. In this respect, we build initial conditions within a more realistic 
technique [HI [12]. Namely, a stationary state with the proper initial conditions is 
produced in an asymmetric trap and then the trap is non-adiabatically transformed to 
the symmetric form. 

The paper is organized as follows. The calculation scheme is sketched in Sec. [2] 
Results of the calculations are discussed in Sec. [3j A summary is given in Sec. [U 



2. Calculation scheme 

We solve the 3D time-dependent Gross-Pitaevskii equation [21] 

ih— (r,t) = [-—V 2 + V(r) +g \*(T,t)\*\*(T,t) (1) 

for the total order parameter \l/(r, t) describing the BEC in both left and right wells of 
the trap. Here go = 4:TTh 2 a s /m is the interaction parameter, a s is the scattering length, 
and m is the atomic mass. The trap potential 

V(r) = ^{ujy + ujy + u; 2 V) + V cos 2 (vrx/g ) (2) 
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includes the anisotropic harmonic confinement and the barrier in x-direction, where Vo is 
the barrier height and go determines the barrier width. Note that the barrier parameters 
may depend on time at the stage of preparation of the initial conditions. 

Following the experiment [HJ[T2], we use a BEC of 87 Rb atoms with a s = 5.75 nm. 
The trap frequencies are u x = 2tt x 78 Hz, uj y = In x 66 Hz, u z = 2ir x 90 Hz, i.e. 
uj y + uj z = 2uj x jTH IT21 . The barrier parameters are V = 420 x h Hz and q = 5.2 /im 
[TTjlT2"]. The distance between centers of the left and right wells is then d =4.4 /xm. For 
iV=1000 atoms, we reproduce the conditions of the JO/MQST Heidelberg experiment 
[TT| [12] for a weak coupling of the left and right BEC fractions through the barrier. 

The static solutions of (OQ) are found within the damped gradient method [27] 
while the time evolution is computed within the time-splitting (28] and fast Fourier- 
transformation techniques. The order parameter ^(r, t) is determined in the 3D 
cartesian grid. The requirement J ' ™ dr 3 \^/(r, t)\ 2 = N is directly fulfilled by using 
an explicit unitary propagator. Reflecting boundary conditions are used throughout, 
but have no impact on the dynamics because of the harmonic confinement. In both 
static and time- dependent cases, only the order parameter with the lowest energy 
(chemical potential /i) is explicitly computed. No time-space factorization is used. The 
conservation of the total energy E and number of atoms N is controlled. Note that the 
Gross-Pitaevski equation is mathematically equivalent to the non-linear Schrodinger 
equation. In this sense, our approach is a counterpart of the time-dependent Hartree- 
Fock method for the system of interacting bosons. 

The JO and MQST are studied in terms of the time-dependent normalized 
population imbalance z and phase difference 9, 

z(t) = NL{t) ~ NR{t) , 6(t) = Mt)-Mt), (3) 

where Ni t R(t) are respectively populations of the left and right wells (with Ni(t) + 
Nn(t) = N) and 4>L,n(t) are the corresponding BEC phases. The populations read 

/+oo 
dr 3 |^(r,i)| 2 (4) 
-oo 

with j = L,R and = *(x < 0,y,z,t), ^ R (r,t) = ty(x > 0,y,z,t). 

The phases 4>j(t) are defined as 

(t>j(t) = arctan (5) 

3 (At) 1 ; 

with the averages 

i r+oo 

7i(«) = W / rfr 3 Im(^.(r,t))|^(r,t)| 2 , (6) 

1 p+oo 

<,■(*) = tt / ^Re^.^t))!^^,^! 2 . (7) 

^ j J -oo 

Computation of the phase time evolution through arctan may be cumbersome. So we 
use (jSJ) only for the static case while the time evolution is calculated through the phase 
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Table 1. Parameters of the procedure for the JO/MQST initial conditions. See text 
for detail. 



N 


d, fim 


r, ms 


q, fim 


JO 


MQST 


JO 


MQST 


JO 


MQST 


1000 


0.25 


0.5 


4 


8 


5.2 


5.2 


3000 


0.45 


1.05 


5 


9 


5.2 


5.2 


5000 


0.62 


1.35 


5 


5 


5.2 


6.6 


10000 


1.35 


1.35 


8 


2 


5.2 


10.1 



increments 4>j(t + St) ~ <f)j{t) + S<fij(t) for a small time step St. Namely, we use 



7?(* + St) + Q{t + St) 



(8) 



with S 1:i (t) = 7j -(t + St) - jj{t), SQ{t) = Q(t + St) - Cj(t). 

The calculations are performed for iV=1000, 3000, 5000, and 10000 atoms. For all 
the cases, the initial population imbalance zq is 0.3 for JO and 0.6 for MQST [TTJ [T2] . 
The initial phase difference is 8^=0 for both JO and MQST. 

In most of the previous studies, the initial state is prepared as the lowest non- 
stationary state with the constrained z in the symmetric trap, see, e.g. [8J. The 
constraint is reasonable for a weak coupling (iV=1000 in our study) but may be 
questionable for SC where the chemical potential [i x associated to the motion in 
x-direction approaches or even exceeds the barrier height. Besides, the constraint 
procedure deviates from the actual experimental initialization [TTJ, [12] of the JO/MQST 
dynamics. Hence we use here, in addition to the constraint calculations, a more 
reliable and realistic initialization following [TTJ [12]. We start from the asymmetric 
trap produced from the symmetric one by a right-shift d of the barrier. The value of the 
shift is adjusted to provide the given z in the lowest stationary state of the asymmetric 
trap. In some cases, an additional widening of the barrier, q — > q, in the asymmetric 
trap is applied. Then the trap is rapidly (for a time r) returned to the symmetric form 
and the JO/MQST time evolution starts. 

The parameters of the procedure are given in Table 1. The return time r is chosen 
to provide a reasonable initialization of JO/MQST. The calculations show that a too 
rapid return shakes the system and leads to a fussy and fragile JO/MQST. Instead, for 
too slowe return, the equilibration process noticeably modifies z and 9 from their initial 
values zq and 6 . Altogether, this procedure is more justified and realistic than the mere 
constraint. For iV=1000, it fully reproduces formation of the initial conditions in the 
experiment [TTJ [12] . 
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3. Results and discussion 

3.1. Static interaction effect 

In our study, we swap from weak to strong coupling (SC) and approach the crossover 
point by increasing the number of atoms N from 1000 to 10000. This results in rising 
the integral interaction effect and a subsequent growth of the chemical potential [l x . 
Thus we naturally come from deeply sub-barrier to above barrier cases. 




0-1 — ■ . , . : — 10 0-1— . , . HO 

-5 5 -5 5 

x [urn] x [urn] 



Figure 1. The double- well trap potential V(x) (bold curve), the linear (without 
interaction) chemical potential fiQ X (dotted line), the non- linear (with interaction) 
chemical potential fj, x (solid line), and BEC density p(x) (dotted curve) in the stationary 
states with N 1 = N 2 = N/2 for 7V=1000 (a), 3000 (b), 5000 (c), and 10000 (d). 

This is demonstrated in Fig. 1 where the static BEC density p(x) and chemical 
potential \i x (for the motion in x-direction) are exhibited at different values of iV 
and compared to the relevant trap potential V x (x) = muj\x 1 jl + Vo cos 2 (7rx/go)- As 
mentioned above, the trap parameters are the same as in the experiment [TH [12]. The 
density 

/+oo 
dydz\^(x,y,z)\ 2 (9) 
-oo 

is determined from the solutions of ([1]) for interacting BEC. 

Figure 1 shows that, for iV=1000, the overlap of the left and right BECs is small. 
The ratio of the densities at their maxima (x = ±2.2 fim) and center of the trap (x = 0) 
is p m / pc—5.1. This is the case of the weak coupling used in the experiment [HJ IT2] . 
Increasing N results in rising the densities and larger overlap in the barrier region. At 
iV=10000, we already have the SC case with p m /p c =2. Here the overlap of the left 
and right BECs cannot be neglected and the system should be treated with the total 
order parameter. At the same time, the left and right density bumps are still well 
distinctive. So a physical view of the system in terms of two (strongly coupled) BECs is 
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still reasonable and one may yet expect for the JO/MQST dynamics described by the 
relative variables z and 9. 

To discriminate the sub-barrier and above-barrier cases, one should compare the 
chemical potential \i x and the barrier height Vo- For fi x , only the motion in x-direction 
has to be taken into account. In the linear case {go=0), we straightforwardly get 

h 

Since uj y +uj z = 2ui x [HI [12], the subtractive term in ( FlOl) is equal to hui x . The calculations 
give a = Uxo/no = 3/4, i.e. just x-motion mainly contributes to the total chemical 
potential fio- This is because the barrier separates the harmonic x- confinement into two 
more narrow regions and thus effectively increases u x . 

In the nonlinear case (go ^ 0), the estimation of fi x is straightforward for ID system 
but complicated for 3D one. Here we roughly put \i x = an where \i is the total nonlinear 
chemical potential. Hence we suppose that in the linear and nonlinear cases the relative 
contribution of x-motion to the chemical potential is the same. 

Figure 1 compares the linear /j, x0 and nonlinear fix to the barrier height Vq. It is 
seen that fixo does not depend on N. It is always much lower than Vo, thus leading to 
the deeply sub-barrier case. If the repulsive interaction is switched on, the nonlinear 
chemical potential \i x rises with N. We see the sub-barrier case for N = 1000, the 
crossover region for iV=3000-5000, and the above-barrier case for iV=10000. 



3.2. JO and MQST evolution 

The evolution of JO and MQST with N is demonstrated in Figs. 2-3. In all cases, the 
initial (t=0) conditions are z =0.3, 9 =0 for JO and z =0.6, 9 =0 for MQST [TT1[T2]. 
Both the constraint technique (CT) and barrier-shift technique (BST) [HI [12], described 
in Sec. [2, are used for initialization of the dynamics. For the BST, the time r when the 
asymmetric trap is fully reduced to the symmetric form is marked by a vertical line. It 
is seen that CT at t=0 and BST at t — r give somewhat different z and 9. The larger N 
and t (and thus the equilibration time), the more the difference. Nevertheless, the CT 
and BST usually initiate a similar (up to a constant time shift) dynamics, especially for 
JO. In what follows, we will mainly exam the BST results. 

First of all, note that for iV=1000 our calculations well reproduce the JO and MQST 
experimental data [HI [12]. Following Fig. 2 a),e), we obtain for JO the robust z- and 
^-oscillations with the frequency ojjo = 2nx23 Hz which is close to the experimental 
value cJjq P = 27rx25 Hz. In Fig. 3a) for MQST, the calculations give ^-oscillations 
around (z)=0.6 with the frequency wmqst = 27rx72 Hz close to the experimental value 
w mqst = 27rx78 Hz. A small underestimation of the experimental frequencies can be 
caused by using in our calculations a smaller number of atoms, iV=1000 instead of 
iV w 1150 ± 150 in the experiment. In Fig. 3e) for MQST, the computed 9 linearly rises 
with the rate 9 = 2ttx 75 Hz^ c^mqst similar to the experimental rate 27rx78 Hz. The 
calculated oscillation amplitudes Az = z max — ^ min =0.6, A9 = 9 max — 9 min =lAiT for JO 
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Figure 2. Time evolution of JO population imbalance z (upper plots) and phase 
difference 6 (bottom plots) for N=1000, 3000, 5000, and 10000, as indicated. The 
barrier-shift (solid line) and constraint (dotted line) techniques are used to initiate the 
evolution. For the barrier-shift case, the return time r is marked by the vertical dash 
line. See text for more detail. 



and Az=0.3 for MQST also reproduce the experimental data. A good agreement of our 
results with the experiment [HI [12] proves high accuracy and realistic character of our 
method and justifies its application to the more sophisticated cases considered below. 

The evolution of JO with iV is exhibited in Fig. 2. It is seen that, despite a 
significant change in the conditions from iV=i000 (weak coupling, sub-barrier transfer) 
through iV=3000-5000 (significant coupling, crossover region) to N= 10000 (strong 
coupling, above-barrier transfer), the JO keep the main features: oscillations of z and 
9 with the same frequency wjo around zero average values. The time shift between z 
and 9 is about one-half a period. Amplitudes of the oscillations do not change with N. 
The evolution with iV (or similarly with the interaction U ~ go between BEC atoms) 
is mainly reduced to a growth of the frequency ojjq. This trend is natural since the 
larger N, the higher the chemical potential /j, x and the larger the barrier penetrability 
K. Furthermore, the larger K, the higher ujo, as it should be in Rabi-like oscillations. 
Altogether we see that JO survive (with a higher frequency) even at SC and above- 
barrier transfer. The crossover region iV=3000-5000 is passed monotonically. 

Note that the above JO evolution is also supported by the weak-coupling arguments 
[9] though application of these arguments needs a word of caution. In the TMA 
weak-coupling picture for the interacting BEC [!9J, JO dynamics is driven by the 
interaction/coupling ratio A = NU/K. With increasing N, the interaction part NU 
grows linearly while the barrier penetrability K rises exponentially. Altogether, A falls 
with which, following TMA calculations [9J, should lead to decreasing cj jo . Instead, 
our calculations demonstrate the opposite trend. The point is that our approach is a 
counterpart of the time-dependent Hartree-Fock method for boson systems, where the 
many-body problem for interacting bosons is reduced to one-body problem for a motion 
of a boson in an effective one-body potential involving the impact of the interaction. 
So the weak coupling arguments [9] should be used in the noninteracting limit where 



Strong-coupling dynamics of Bose-Einstein condensate in a double-well trap 



9 



cjjo oc K. Hence the JO trend in Fig. 2. 

In Figure 3, the evolution of MQST with N is demonstrated. As compared to JO, 
this evolution is more complicated and needs more time to be exhibited. Hence we use 
here the larger time interval 80 ms. Fig. 3 shows that MQST is transformed with N 
to JO from Fig. 2. Namely, at iV=1000, there is an ordinary MQST in agreement with 
the experiment [HJ[12|. At iV=3000-5000, z-oscillations around (z)=0.6 are gradually 
reduced to slower oscillations around (z)=0 and, at N= 10000, basically converge to JO 
in Fig. 2d). The linear evolution of 9 is turned into JO-like oscillations around the 
average {6)=2tt. Since the shift 2tt is irrelevant, one actually gets the JO in Fig. 2h). 




Figure 3. The same as in Fig. 2 but for MQST dynamics. 



A similar MQST evolution was earlier obtained in the strong interaction limit for ID 
BEC [26] . It was shown that at sufficiently strong interaction the amplitude of MQST z- 
oscillations starts to grow with NUm, thus leading to reappearance of tunneling between 
two wells. This effect was explained, within the modified TMA, by coupling the second 
and third modes. Our model takes into account all the modes and thus should cover 
the result [26] but now for 3D system. Actually, the transfer MQST — > JO in our 
calculations exhibits the similar effect: increasing the amplitude of ^-oscillations and 
thus reappearance of the tunneling. 

The most remarkable result of the present study is that, despite a significant overlap 
and strong coupling in the SC case, the left and right BECs in the trap still keep 
their individuality and accept the Josephson-like dynamics in terms of the relative 
variables, population imbalance z and phase difference 9. In particular, the JO dynamics 
obtained at a modest initial z Q remains very regular and robust. This means that 
BEC interferometry and related processes may be successfully realized not only at weak 
coupling with a deeply sub-barrier transfer but also at SC with above-barrier transfer. 
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4. Summary 

The dynamical evolution of coupled Bose-Einstein condensates (BEC) in a double well 
trap was investigated while modifying the system from a weak coupling case (small 
overlap of BECs, deeply sub-barrier transfer) to a strong coupling case (considerable 
overlap of BECs, above-barrier transfer). The evolution was driven by increasing the 
number N of BEC atoms and thus rising the total effect of the interaction between 
BEC atoms. The numerical analysis was performed by solving the 3-dimensional time- 
dependent Gross-Pitaevskii equation. Thus the two-mode and many other ordinary 
approximations were avoided. The main dynamical regimes, the Josephson oscillations 
(JO) and Macroscopic Quantum Self- Trapping (MQST), were inspected. 

The calculations show that the JO successfully survive even at strong coupling but 
acquire a higher frequency. The MQST is destroyed and finally reduced to the same, 
though more fragile, high-frequency JO mode. Altogether we see that the Josephson- 
like dynamics certainly persists and remains robust at strong coupling. This means 
that, despite a strong overlap, the left and right BECs in the trap still keep their 
individuality and the relative variables, population imbalance z and phase difference 9, 
remain reliable. These findings show that BEC interferometry and related phenomena 
may be extended to the case of strong tunneling coupling, which opens a new avenue 
for further explorations. 
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